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ABSTRACT 

A significant fraction of unstable multiple planet systems are likely to scatter during 
the transitional disc phase as gas damping becomes ineffectual. Using a large en- 
semble of FARGO hydrodynamic simulations and MERCURY n-body integrations, we 
directly follow the dynamics of planet-disc and planet-planet interactions through the 
clearing phase and on through 50 Myr of planetary system evolution. Disc clearing 
is assumed to occur as a result of X-ray driven photoevaporation. We find that the 
hydrodynamic evolution of individual scattering systems is complex, and can involve 
phases in which massive planets orbit within eccentric gaps, or accrete directly from 
the disc without a gap. Comparing the results to a reference gas-free model, we find 
that the n-body dynamics and hydrodynamics of scattering into one- and two-planet 
final states are almost identical. The eccentricity distributions in these channels are 
almost unaltered by the presence of gas. The hydrodynamic simulations, however, also 
form a population of low eccentricity three-planet systems in long-term stable config- 
urations, which are not found in n-body runs. The admixture of these systems results 
in modestly lower eccentricities in hydrodynamic as opposed to gas- free simulations. 
The precise incidence of these three-planet systems is likely a function of the initial 
conditions; different planet setups (number or spacing) may change the quantitative 
character of this result. We analyze the properties of surviving multiple planet systems, 
and show that only a small fraction (a few percent) enter mean-motion resonances af- 
ter scattering, while a larger fraction form stable resonant chains and avoid scattering 
entirely. Our results remain consistent with the hypothesis that exoplanet eccentric- 
ity results from scattering, though the detailed agreement between observations and 
gas-free simulation results is likely coincidental. We discuss the prospects for further 
tests of scattering models by observing planets or non-axisymmetric gas structure in 
transitional discs. 

Key words: planets and satellites: dynamical evolution and stability — planet- 
disc interactions — planetary systems — protoplanetary discs — hydrodynamics — 
scattering 



1 INTRODUCTION 

Prior to gap opening, gravitational interactions between 
planets and their surrounding gas discs act to efficienctly 
damp planetary eccentricity. This damping implies that the 
formation of giant planets via core accretion leads to initially 
circular orbits, and that post-formation dynamical effects 
must be sought to explain the observed eccentricity distri- 



bution of massive extrasolar planets (Butler et al 



Dynamical instability in multi-planet systems, leading to 
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orbit crossing and the ejection or merger of some of the 
planets, is one way to achieve this (Rasio & Ford 1996] [Wei^ 
[denschilling Marzari|l996| pn &; Ida|1997| ). Extensive n- 
body experiments have shown that with realistic mass func- 



tions, the dynamics of unstable two-planet (Ford & Rasio 



2008), three-planet (jChatterjee et al. 2008), or richer sys- 



tems (Juric & Tremaine 2008) can successfully reproduce 
the observed eccentricity distribution. 

The good agreement between n-body dynamical experi- 
ments and observations is somewhat surprising, as it is hard 
to envisage a consistent scattering scenario in which gas disc 
interactions are not important. The simplest model is one 
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in which the typical planetary system forms with just two 
massive planets. In this case, the only unstable systems have 
planets that are separated by less than the threshold for 

11993) 



Hill stability (Marchal & Bozis 



1982 



Gladman 



and 

most of the unstable subset scatter on a time scale that is 

n5 



much shorter than the gas disc dispersal time (~ 10^ yr 
Wolk Walter|1996 ). Scattering is therefore likely to occur 
as soon as gas disc damping declines below some thresh- 
old level, and will be synchronized with the epoch of disc 
dispersal. Related arguments apply, less robustly, to richer 
planetary systems. For more than two planets there is no 
analytic criterion for stability, and relatively loosely packed 
initial configurations can delay the onset of scattering for 
tens or hundreds of Myr (Chatterjee et al.,,2QQ8J . Absent 
fine tuning, however, it is not possible to construct initial 
conditions in which all of the instabilities occur on > Myr 
time scales, without over-producing stable low eccentricity 
systems that are not seen observationally. As the gas dis- 
sipates there should therefore be a significant number of 
systems primed to go unstable, scattering in the presence of 
at least some remnant gas. 



The large computational cost of hydrodynamics rela- 
tive to gravity has necessarily motivated approximate ap- 
proaches that marry disc physics to n-body simulations 
of scattering. One method is to include an analytic de- 
scription of planet-disc interactions as a force within an n- 
body calculation. This has been used to study the migration 
of planets into resonance (e.g. Lee &: Peale|2002" Snellgrove 
et al.||2001 ), and through resonance into scattering (Adams 
fc Laughliii|2Q03||Mo orhead & Adams|2005l|Lee et al.|2009| ). 
Scattering with more than two planets in an embedding disc 



has been studied by Chatterjee et al. (2008) and Matsumura 
et al. (2010), with the disc modeled in ID and torques on 



the planets calculated from that density profile. There have 
also been preliminary studies that include a full hydrody- 
namic treatment of the gas disc during active periods of 
planet scattering. Moeckel et al. (2008|) performed 3D sim- 



ulations starting from very unstable equal-mass, two-planet 
systems with a remnant gas disc exterior to the planets' or- 



bits. Marzari et al. (2010) performed 2D simulations study- 



ing migration and resonant trapping in three-planet systems, 
some of which led naturally to scattering when resonance 
was broken. While both of these studies hinted at possibly 
significant dynamical effects mediated by the gas, they were 
too limited in scope to attempt statistical comparison to 
gas-free results. 

Our goal in this work is to directly quantify the impact 
of residual gas discs on the outcome of planet scattering. 
We assume that the interesting dynamics is coincident with 
the epoch of disc dispersal, and follow the planet-planet 
and planet-disc interactions through this transition phase 
using a 2D hydrodynamic code. We consider three-planet 
systems, setup almost identically to previous gas-free calcu- 
lations d Chatterjee et al.|[20Q8l [Raymond et al.||2010| ), and 
model disc dispersal using a physically motivated photoe- 
vaporation model (Owen et al. 2011). The hydrodynamics 
are run long enough to take the systems through the tran- 
sition phase into the pure n-body regime, after which we 
switch to a standard integrator for a further 50 Myr of grav- 
itational interactions. 



2 THE SIMULATIONS 
2.1 Initial conditions 

To facilitate direct comparison with previous n-body results, 
we consider three-planet systems around stars with mass 
= 1 Mq, with the mass m of each planet drawn ran- 
domly from a power law distribution f{m) oc m~^'^ in the 
range 0.3-5 Mj. The innermost planet is placed with a semi- 
major axis ai = 3 au, and the spacing of the other plan- 
ets is determined following the approach of [Chambers et ah] 
(1996), with each adjacent planetary pair separated by a 
fixed number of mutual Hill radii Rnij- That is, with 



RHij = 



rrii + rrij 
3M. 



1/3 



ai + aj 



(1) 



the middle planet is placed at a2 = ai + KRh 12, and the 
outer planet at as = a2 + KRh22,- The constant K largely 
determines the timescale on which gravitational instability 



2002 



occurs (Chambers et al. 1996 Marzari & Weidenschilling 



Chatterjee et al. 2008). We set K at approximately 



the maximum value that does not overproduce stable sys- 
tems, whose low final eccentricity would be inconsistent with 
exoplanet observations. With our mass spectrum, experi- 
mentation showed that = 4.0 yielded a distribution of 
instability times with a median ^ 4 x 10^ yr, and ~ 10% 
of the setups unscattered after 50 Myi[^ 

The planets are initially embedded in a gas disc which 
we model between 1 and 40 au. The surface density initially 
follows a radial profile 5](r) = 200 (r/1 au)~^ g cm~^. We 
include photoevaporation of the disc using the mass-loss pro- 
files of Owen et al. (2011), which are based on the coupled 



hydrodynamic and radiative transfer models introduced by 



Owen et al. (2010). In these models. X-ray luminosity drives 



photoevaporative disc dispersal. Our discs initially contain 
just over 5 Mj of gas, which is both the typical total mass 
of the three planets and the mass at which the photoevap- 
orative mass loss begins to alter the disc structure on a 
fast timescale. We pragmatically set the X-ray luminosity 
of our star at 1.4 x 10^° erg s~^, which yields a disc clearing 
timescale of just over 10^ years depending on the details of 
the planet-disc interactions (see Figure [T]). This timescale, 
coupled with our chosen planet spacing, means that about 
half our systems might be expected to scatter when there is 
still a significant amount of disc material left. 

We evolve 100 random realizations of the planetary ini- 
tial conditions within the fixed initial conditions for the disc. 
We use direct hydrodynamic simulations to model the cou- 
pled planetary and disc dynamics while the disc mass is high 
enough to still influence the planets' evolution, followed by 
a longer n-body run once the disc has been cleared. For 
the hydrodynamic portion of each run we use a modified 



version of the 2D code FARGO ( |Masset[[2000|, and for the 
-body foUow-on we use MERCURY (| Chambers|1999| . 

FARGO is a 2D Eulerian hydrodynamic code tailored to 



^ This instability timescale is about an order of magnitude larger 
than the result for K = 4.0 with equal mass planets found by 
[Chatterjee et al.| ( [2008| >, where for a given K each planetary pair 
is at a fixed separation relative to the major mean motion reso- 
nances. With our mass spectrum and K = each adjacent pair 
is somewhere between roughly the 3:2 and 2:1 mean motion res- 
onances. 
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Figure 1. The evolution of the remaining disc mass for each 
run is plotted in gray. The red histogram shows the (arbitrarily 
normalized) distribution of dissipation times, i.e. when the disc 
mass reaches 0. 



planet-disc interactions (Masset 2000). We have modified 



the code somewhat to make it more suitable for this work 



(similar modifications were made by Marzari et al. 2010). 



The n-body integrator was changed from a 5th order Runge- 
Kutta to a (7,8) order embedded pair (Prince & Dormand 
1981). We include stepsize control in the following way: the 



results of the 7th and 8th order estimate of each compo- 
nent of each planet's velocity and position are compared at 
the end of each hydrodynamic timestep. If the fractional 
difference between the two exceeds 10~^, the planets are 
integrated through the hydrodynamic step with a smaller 
n-body timestep, repeating this process until our accuracy 
criterion is met. 

We have also implemented the removal of escaping plan- 
ets if they are on a hyperbolic orbit and more than 400 au 
from the star, as well as collisions between planets and be- 
tween planets and the star. Planet-planet mergers are mass 
and momentum conserving and take place whenever two 
planets' surfaces are detected to touch, with all planets as- 
sumed to have radii of 1.3 Rj. Planets merge with the star 
if they approach within 5 Rq. 



2.2 Numerical details 

Photoevaporation is included using the X-ray driven mass- 
loss profile of 



Owen et al. 



(2011), scaled to our 1 Mq star 



and chosen X-ray luminosity of 1.4 x 10*^^ erg s~^. The az- 
imuthally constant but radially dependent mass loss rate 
from these profiles is imposed across the numerical grid as 
an externally specified sink term in the continuity equation. 
In these models photoevaporation opens a gap near 2 au, 
and when the gas interior to the gap is gone the directly- 
illuminated outer disc erodes from the inside out (see Figure 
9 of Owen et al.||2011 ). We implement this transition to di- 
rect illumination somewhat approximately, as we do not fol- 
low the disc inwards of 1 au and set a lower surface density 
limit of 2 X 10~^ g cm~^ below which we do not photoevapo- 
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Figure 2. Resolution tests in which a 1 Mj planet is placed on 
an eccentric orbit in an initially unperturbed disc. The series of 
blue lines shows a disc with 11(1 au) = 200 g cm~^, and the red 
lines have E(l au) = 50 g cm~^. In each series the lightest shade 
is run with ne = 256, the medium shade with no = 512, and the 
dark shade with ng = 1024. 



rate any material for numerical reason^ As the simulation 
progresses, we track the azimuthally averaged surface den- 
sity of each ring. When the innermost radius has an average 
surface density below 8 x 10"*^ g cm~^, we rescale the inner 
radius of the photoevaporation profile to that radius. 

We allow accretion onto the planets using a similar ap- 
proach to the unmodified version of FARGO, based on the 



prescription used by Kley (1999). The planet's Roche lobe 



(approximated by a constant radius Rr) is drained on some 
timescale, with material inside 0.45 Rr drained at twice the 
rate of that between 0.45 and 0.75 Rr. Material outside 0.75 
Rr is left alone. The timescale on which the Roche lobe is 
drained is set to twice the planet's orbital period. With our 
chosen parameters, the disc is typically dispersed via a com- 
bination of planetary dynamics and photoevaporation (with 
a small contribution from accretion) between 1-1.5x10^ yr, 
with a few instances stretching to 0.8-1.9x10^ yr. In Figure 
[T]we show the mass evolution of all the discs in our runs. 

^ With this density floor, our flnal disc mass when photoevapo- 
ration has run its course is ~ 10""^ M /. 
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Figure 3. The evolution of the planetary system and the gas disc for an example run. The orbits of the planets are plotted on top of 
the surface density plots in the same colours as the plots of the orbital elements and masses. The surface density plot on the left shows 
the system before scattering, when the red, green, and blue planets have partially cleared a mutual gap. During the scattering event the 
green and blue planets collide to form the orange planet on an orbit with e ~ 0.6. The aftermath of the scattering is shown in the middle 
density plot. Over the next ~ 2.5 x 10^ years, disc interactions damp the eccentricity of the outer planet; after disc clearing, the two 
planets are stable. The right density plot shows the system near the end of the photoevaporative clearing, with each planet orbiting in 
a gap. 



We use open boundaries at the inner and outer edges 
of the grid. We use a fixed disc aspect ratio of H/r — 0.04, 
and an alpha-viscosity with a — 10~^. The grid is setup 
with ne — 512 azimuthal cells, and Ur — 207 rings. The ra- 
dial spacing is logarithmic, with approximately square cells 
between 2.5 and 10 au, but coarser radial spacing between 
1-2.5 and 10-40 au so that those cells are rectangular with a 
2:1 aspect ratio. Our highest resolution is thus concentrated 
in the region of highest interest, where the majority of the 
planets spend most of the time. 

We tested the adequacy of the grid resolution by setting 
up a planet with m = 1 Mj, and eccentricity e = 0.5, in 
an initially unperturbed disc. We considered two initial disc 



masses, S(l au) = 200 and 50 g cm~^, and three resolutions, 
no — 256, 512, and 1024. The results are shown in Figure 
|2] The more massive disc results are shown in blue, and the 
low mass disc in red. The lightest shade is ne — 256, and 
the darkest shade is ne — 1024. The low mass disc is well 
converged in both the orbital evolution and planet mass. 
The high mass case is only moderately converged. While the 
orbit of the planet ends up at approximately the same place 
for all three resolutions, the mass evolution of the planet is 
quite different. The differences in the orbital evolution are 
most pronounced at the time when the planet is accreting 
the most mass, between 1000-2000 yr. Visual inspection of 
the surface density of these runs shows that gas is accreting 
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Figure 4. The evolution of the planetary system and the gas disc for an example run. The orbits of the planets are plotted on top of the 
surface density plots in the same colours as the plots of the orbital elements and masses. The surface density plot on the left shows the 
system before scattering, when the red, green, and blue planets have partially cleared a mutual gap. The middle density plot is during a 
period of active planetary dynamics. After this stage the disc damps the eccentricity of the three planets, which remain in roughly this 
configuration through 10^ yr. This phase is shown in the right density plot. 



onto the planet in streams that are not well resolved at low 
ne. We conclude that for studying the orbital evolution of 
the planets — the prime focus of this work — ne = 512 
is sufficient. Higher resolution, along with a more physical 
treatment of mass accretion within the Roche lobe, would 
likely be needed to model mass evolution in detail. 

The results of the resolution test suggest that a massive 
planet scattering into an undisturbed disc will accrete ma- 
terial at a resolution-dependent rate until a gap is opened. 
This manifests itself as an offset in the mass evolution of the 
disc and the planet (see the similar slopes of the Mpianet and 
Mdisc curves in Figure |2] after ^ 1500 years). Runs that scat- 
ter very early, when a large amount of gas is still in the disc, 
are susceptible to this numerical effect, which may alter the 
disc dissipation timescale. Offsets of the order 0.1 Mj will 



not move the disc mass evolution out of the main band seen 
in Figure [l] though, and are certainly small compared to the 
expected physical dispersion in dissipation history that re- 
sults from varying stellar X-ray luminosity. Hence, we do 
not consider this to be of great importance. 

Each set of initial conditions is integrated with FARGO 
until the disc is photoevaporated down to our density floor, 
at which point the run is terminated. The planetary system 
is then evolved further using the FARGO integrator with no 
disc until 2 x 10^ yr. At this point, the system is turned 
over to the Burlisch-Stoer integrator of mercury (Cham- 
1999| and integrated to 50 Myr. As in the FARGO 



bers 



portion of the integration, planetary mergers and ejections 
are allowed with the same parameters. These hydrodynamic 
runs are then compared to an ensemble of 500 purely n- 



© 2011 RAS, MNRAS 000,[T}{T2] 



6 A^. Moeckel and P. J. Armitage 



^00[ 




2 (not stable) 

nD=3 



nn=2 



nn=1 



1 2 3 4 5 6 7 8 

log (ti / yr) 



disc no disc 



Figure 5. Left: The distribution of the instabihty times for the 
runs with a disc (red) and without (black). The final bins include 
all runs still stable at 50 Myr for the disc runs, and 15 Myr for the 
discless runs. RightThe final state of the planetary systems after 
50 Myr (for the disc runs) and 15 Myr (for the discless runs). 
The bins labeled 'rip = 1' and 'rip = 2' refer to single planet 
systems and double systems that are formally Hill stable. The 
bin labeled '2 (not stable)' refers to 2 planet systems that are not 
Hill stable. Triple planet systems ('rip = 3') include only those 
that underwent a scattering event but retain three planets at the 
end of the simulation. 



body simulations, set up using identical initial conditions 
but integrated using mercury throughout. 



3 RESULTS 

The hydrodynamic runs are characterised by several phases: 
an initial clearing of gaps and relaxation of the disc, the 
onset (or not) of dynamical instability and planet scattering, 
renewed relaxation of the disc and eccentricity damping, and 
clearing of the disc by photoevaporation. These phases are 
illustrated for two example runs in Figures [3] and ^ In each 
of these examples scattering occurs, and the once-prominent 
gaps are filled in. As the disc comes to terms with the new 
planetary configuration eccentric gaps develop (also noted 
by Marzari et al.|20To ) , and finally the disc is cleared leading 
to the purely gravitational phase of our experiments. Here 
we discuss the most interesting aspects of these simulations. 



3.1 Instability timescale and outcome 

In Figure |5] we compare the distribution of instability times 
between the hydrodynamic and n-body control runs. Dur- 
ing the FARGO runs, the time of scattering is identified as 
the first time the semi-major axis of any planet changes by 
10 percent between two output times. Visual inspection of 
the results confirms that this is a reasonable estimate of the 
onset of instability. During the mercury runs, an encounter 
between two planets within 3 Hill radii is the criterion. The 
bulk of the distribution is similar for both sets of simulations; 
the main difference is at the long-time tail of the distribu- 
tion, where the FARGO runs have roughly 50% more runs 
that never go unstable over the course of the simulations. 

In the right panel of the figure we plot the percentage of 
unstable runs ending in a given planetary configuration. The 
main difference between the simulation sets is the number of 
unscattered systems and the presence, in the hydrodynamic 
runs only, of a population of scattered systems that retain 
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Figure 6. The eccentricity distribution at 10^ yr for observed 
exoplanets (blue), the runs with no disc (black), and the runs 
with a disc (red). The dashed line shows the disc-less result with 
random inclinations up to 1°. Top panels show only single and 
Hill stable double systems: K-S and A-D tests do not conclusively 
rule out the null hypothesis that the latter two distributions are 
the same (p- values 0.10 and 0.037, respectively). Bottom panels 
include all systems that scattered prior to 50 Myr, i.e. including 
triples and non Hill stable doubles. K-S and A-D tests reject the 
null hypothesis for these samples (p- values < 10~^) 
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Figure 7. Semi-major axis-eccentricity scatter plots of the stable 
systems at 10^ yr with and without the disc. A two-dimensional 
K-S test rejects the null hypothesis that the two distributions are 
the same (p- value 0.041). 



all three planets. There are also a handful of double systems 
that are not provably Hill stable in the mercury runs, but 
not enough to alter any conclusions. We discuss the differ- 
ence in the percentage of unscattered triple systems further 
in sections [3.31 and [3^ 
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Figure 8. The spacings of the triple systems in the disc runs 
(the 'rip = 3' bin in Figure [Sj. Each row shows one system: smah 
grey points show the initial spacing of the planets, and the large 
coloured points show the semi-major axes at 50 Myr. The colours 
correspond to the initial ordering: red points were originally the 
innermost planet, green points were in the middle, and blue points 
were outermost. The numbers between adjacent points show how 
many symmetric Hill radii separate the pair. The number to the 
right of each set shows the scattering time taken from extended 
runs to 10^ yr; only two systems went unstable in that time. 



3.2 Eccentricity distributions 

With realistic planet masses, scattering from nearly (but not 
identically) coplanar and circular initial conditions is quite 
successful at reproducing the observed eccentricity distribu- 
tion of giant extrasolar planets. Good agreement can be ob- 
tained irrespective of the initial number of planets (Ford fc| 
[Rasio|2008i Cha tterjee et al.| 2QQ8; Juric & Tremaine 2008^ 
These authors compare their results to the observed sam- 
pl^ In our study, a subtlety arises from the fact that FARGO 
is a 2D code, and thus we simulate the planetary dynamics 
strictly in the plane. This has some consequences for the ec- 
centricity distribution obtained by purely gravitational stud- 
ies; namely, the distribution is shifted toward lower eccen- 
tricities, clearly seen in Figure [6] We address this point fur- 
ther in appendix [A] When comparing our FARGO runs to a 
disc-less n-body result we thus compare to a coplanar suite 
of experiments run with MERCURY rather than directly to 
the observed sample. With non-zero inclinations, these ini- 
tial conditions yield a distribution quite close to the observed 
exoplanets in a pure n-body suite of experiments, which we 
verify by running the same initial conditions with random 
inclinations up to 1 degree. For the observed sample we ex- 
clude expolanets with a < 0.1 au to take tidal circularisation 
out of consideration. 

To statistically compare the sets of simulations we per- 
form two-sample Kolmogorov-Smirnov (K-S) tests and two- 
sample Anderson-Darling (A-D) tests on the cumulative 
distributions. The latter test is more sensitive to deviations 
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Figure 9. The eccentricity distribution of stable systems (top) 
and all scattered systems (bottom), with the runs with discs di- 
vided into those that scattered with disc material present and 
those that scattered after disc dispersal. K-S and A-D tests re- 
ject the null hypotheses that the early scattering systems follow 
the same distribution as the late scatterers (p-values 0.0016 and 
3.2 X 10-^) or the disc-less runs (p-values 0.0069 and 6.6 x 10"^). 
The late scatterers are ambiguously distinct from the disc-less 
runs (p-values 0.067 and 0.033). Two-dimensional K-S tests on 
the semi-major axis-eccentricity distributions of the early and 
late scatterers likewise reject the null hypothesis only when unsta- 
ble systems are included (p- value 0.071 with only stable systems, 
0.0037 with unstable systems). 



at the tails of the distribution. We reject the null hypothesis 
that the two samples are drawn from the same distribution 
when p-values are less than 0.05. In Figure [6] we show the 
eccentricity distributions of the disc and disc-less cases. The 
top panels show the results including only those systems 
that are provably stable, i.e. single planet and Hill-stable 
double planet systems. The null hypothesis is not ruled out 
by the K-S (p- value = 0.10), though the A-D test rejects it 
(p- value = 0.037). The greater weight the A-D test gives to 
the distributions' tails validates the visual impression that 
the main difference is at the high eccentricity end. For com- 
parison we also plot the observed eccentricity distribution, 
as well as the result of our simulations with inclinations up 
to 1 degree. While the eccentricity distributions alone are 
quite similar for stable systems, the combined eccentricity- 
semi-major axis distribution is marginally different. Figure 
|7| shows scatter plots of these distributions. A 2D K-S test 
(using the method in Fasano & Franceschini 1987) rejects 
the null hypothesis, with a p- value of 0.041. 

The bottom panel of Figure [5] shows the same compar- 
ison between the hydrodynamic and n-body integrations, 
but now including all systems that scattered prior to 50 
Myr, regardless of provable stability. This means that for 
the FARGO runs we include those triples that scattered and 
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Figure 10. The evolution of the planetary system and the gas disc for an example run. The orbits of the planets are plotted on top of 
the surface density plots in the same colours as the plots of the orbital elements and masses. The surface density plot on the left shows 
the system before scattering, when the red, green, and blue planets have partially cleared a mutual gap. During the scattering event the 
green and blue planets collide to form the yellow planet on an orbit with e ~ 0.6. The aftermath of the scattering is shown in the middle 
density plot. Over the next ~ 2.5 x 10^ years, disc interactions damp the eccentricity of the outer planet; after disc clearing, the two 
planets are stable. The right density plot shows the system near the end of the photoevaporative clearing, with each planet orbiting in 
a gap. 



re-circularised in more widely spaced configurations (such as 
the one in Figure [4|, while for the mercury runs we include 
the few double systems that are not provably stable. These 
distributions are visually more distinct, and this is born out 
by the statistical tests, with p- values for both tests less than 
10"*^. To test that the scattered triples are stable on longer 
timescales, we integrated the 12 cases furtheiQ to 10 yr. 
The spacing of these systems at 50 Myr, and the outcome 
of the further integration, are shown in Figure [S] With two 
exceptions, the systems remain stable for 10^ yr. 

Since the disc is the driver of any difference that might 
arise between the distributions, we split the sample into 
those in which instability begins when disc material is still 
present, and later scatterers. These distributions are shown 
in Figure [9] Whether or not we include the unstable sys- 
tems, the distribution of early and late scatterers are signif- 
icantly different, with K-S and A-D p- values of 0.0016 and 
3.2 X 10~^. While the early scatterers are likewise convinc- 
ingly statistically distinct from the disc-less runs (p-values 
0.0069 and 6.6 x 10~^), the two tests also ambiguously reject 
the null hypothesis for the disc-less runs and the late scat- 
terers (p-values 0.067 and 0.033). We would not necessarily 

For these integrations we used the Hybrid integrator of MER- 
CURY, with a 20 day stepsize. 



expect any difference here, as there is no disc to alter the dy- 
namics, and Chatterjee et al.| ( |200 8) showed that the time of 
instability is unimportant in determining the final outcome 
of a scattering experiment. Note that we suffer somewhat 
from low numbers of runs in this late scattering category, 
due to the large number of runs that never scatter over the 
length of our integrations. 



3.3 Mean Motion Resonances: 2 planets 

Convergent migration in a disc can lead to trapping into 
(e.g 



resonance 



Lee Peale 2002; Snellgrove et al. 2001 



Kley et al. 2004 Moor head & Adams 2005), and is seen 



with high frequency in the previous work most relevant to 
this one, the coupled n-body -ID disc models of |Matsumura| 
|et al.| ( 20101 . We first searched for mean motion resonances 
(MMRs) among the 39 stable double systems by observing 
the apsidal alignment vj2 — vji and the resonance angles for 
{p -\- q) \ q MMR given by (e.g. [ Murray Dermott|1999t 

^1,2 = (p + q)M - pAi - qzAJi^2' (2) 

Here subscript 1 and 2 refer to the inner and outer planet, 
A is the mean longitude, and m the longitude of pericenter. 
We searched for resonances up to 8th order by monitoring 
these angles for libration around fixed values. 
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We find just two double systems in 3:1 MMRs witii li- 
bration semi-amplitudes of ^ 150 and ^ 80°. The case with 
the larger libration amplitude did not scatter until the disc 
had been photoevaporated. The low frequency of 2 planet 
resonances, and the large libration angles, are symptomatic 
of random entry into resonance due to purely gravitational 



effects (Raymond et al.|2008), rather than systematic driv- 



ing into resonance due to gas disc migration (Matsumura 
et al.|20T0 ). We observe a dramatically lower resonant frac- 
tion than the ^ 75 percent seen by Matsumura et al. (2010). 



The difference is probably due to the short lifetime of the 
discs in this study compared to the Myr dispersal timescales 



assumed by Matsumura et al.|(2010), which allows less time 



for post-scattering migration to move planets into resonance. 

As noted in section |3.1| the disc runs have a higher 
fraction that remain stable for the entire simulation than 
the n-body only runs. Because our initial conditions natu- 
rally place the planets near to several MMRs, there is the 
possibility that two adjacent planets may get trapped into 
resonance, a situation that can enhance the stability of a 



three planet system (Fabrycky & Murray-Clay 2010). We 
also search each adjacent pair of planets in the systems that 
never scattered for MMRs. Out of 27 such cases, we find 
three in which the inner two planets are in resonance (two 
3:2 and one 2:1), and three cases in which the outer two 
planets are in a 2:1 resonance. These resonant systems are 
potentially part of the reason for the large number of un- 
scattered systems. 



3.4 Mean Motion Resonances: 3 planet resonant 
chains 

Further investigating the unscattered runs, we discovered 
that 12 of the 27 unscattered runs are in three body res- 
onance. In eight cases the inner-middle and middle-outer 
pair move into resonance with librating resonance angles 
(equation [2| , and in three of these the inner-outer pair are 
likewise in a resonance. The chains we find are 9:6:4 (four 
cases), 6:3:2, 3:2:1 (two cases), and 4:2:][^ In each case a 
Laplace angle of the form 



(pj^ — mXi — (m + n)X2 + nAs 



(3) 



is identified (Aksnes 1988), which librates around a fixed 



value for small integer values of m and n, with subscripts 
1,2,3 referring to the inner, middle, and outer planets. 

The other four of the 12 triple resonances are less deeply 
in resonance; all three exhibit libration of the m = 1, n = 3 
Laplace angle but the adjacent planet pairs are not in the 
2:1 resonances. In Figure we plot the resonant behavior 
of three example chains; a representative of the 9:6:4 cases 
(i.e. a double 3:2 resonance), the 6:3:2 case, and the Laplace 
resonance case, the former two of which exhibit asymmet- 
rical resonance in some of the angles. In the 9:6:4 run each 
adjacent pair of planets is trapped into a 3:2 resonance at 
^ 4000 yr, the same time that the Laplace angle begins 
librating around = tt. The inner-outer pair move into 
apsidal alignment, but the resonant angles continue to cir- 
culate. 



The 6:3:2 case moves into triple resonance at ^ 1000 yr, 
roughly the same time as the inner-middle pair lock into a 
3:2 resonance, but somewhat before the other two pairings 
fully enter resonance at ^ 10^ yr. The initial inner-middle 
resonance and the Laplace angle all appear to be settling 
toward resonant angles of or tt, but when middle-outer 
and inner-outer resonance begin, most of the angles move 
into asymmetrical resonances, including the Laplace angle 
which shifts to librate around a value (very) slightly more 
than TT. 

The double 2:1 Laplace resonance shows libration oi (pL 
around an angle completely incommensurate with tt. In this 
case the inner-middle 2:1 pair immediately enters asymmet- 
ric resonance, along with the Laplace angle. The middle- 
outer 2:1 pair are in a symmetric resonance, and the inner- 
outer pair are flirting with asymmetric resonance, though 
continues to circulate. We note that this asymmetric behav- 
ior in the Laplace resonance is a preferred outcome of sim- 



ulations of migration into double resonance ( Beauge et al. 
2008| ), although this particular case did not undergo very 



much migration. It would appear that the high incidence of 
triple resonances, as well as the perhaps enhanced stability 
due to a single adjacent pair migrating into resonance, are 
the underlying reason behind the large number of unscat- 
tered systems in our hydrodynamic simulations compared 
to the n-body runs. 



3.5 Periods of observability 

Disc-planet interactions may be able to excite eccentricity 
to modest values, e 



^ The Laplace resonance, as seen in GJ 876 ([Rivera et al. |2010| ) 
and potentially HD 82943 ( [Beauge et al.|2008| >. 



0.1-0.2 ( [D'Angelo et al.|[2006[ ), but 

planetary dynamics are probably necessary to reach higher 
values. Some of the more interesting features in the disc 
surface density in our simulations occur when eccentricities 
are high, during or just after scattering (see Figures [s] and 
|4|. These include non-axisymmetric features, new gaps be- 
ing opened, and eccentric gaps. If scattering (or at least its 
signature in the eccentricity of the planets) is reasonably 
common in transitional discs, then these features may be 
observable. We estimate the likelihood that these features 
might be observed by calculating the amount of time that 
systems spend with at least one planet's eccentricity above 
a threshold value, while also having at least 1 Mj of gas 
left. This distribution is shown in Figure pT] for three values 
of the threshold eccentricity, emin = 0.1, 0,25, and 0.5. 

For the higher and more interesting thresholds, most 
of the runs do not have eccentric enough planets for any 
reasonable amount of time. If we require that interesting 
features are present for something like 25% of our typical 
disc lifetime, tobs ^ 3 x lO'^ yr, a threshold of emin = 0.25 
means that 10-15% of the runs will have planets in eccen- 
tric enough orbits to be notieceable; thus a few percent of 
transition discs might be expected to harbor planets caught 
in the act or the aftermath of scattering. 



4 CONCLUSIONS 

If planets exit the gas disc phase in non-resonant configu- 
rations, a significant fraction of unstable systems are likely 
to scatter roughly contemporaneously with the transition 
between gas-rich and gas-poor conditions. We have studied 
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Figure 11. The cumulative distribution of tobs 7 defined as the 
amount of time a system with at least 1 M j of gas left in the disc 
harbors at least one planet with e > Cmin- 



whether the presence of small masses of gas, in photoevap- 
orating transition discs, significantly modifies the outcome 
of scattering. Previous work that has addressed this ques- 
tion has been ambiguous. Detailed statistical studies, using 
approximate hydrodynamic treatments of the disc, have sug- 
gested that gas free initial conditions are a reasonable ap- 
proximation to the more physically complete case that in- 
cludes gas ( Chatterjee et al.||2008 Matsumura et al.]|2010 ). 
On the other hand, earlier hydro dynamical work suggested a 



possibly significant effect (Moeckel et al. 2008; Marzari et al, 



2010), but was hampered by very limited statistics. The im- 
mediate aftermath of a scattering event, when multiple grav- 
itationally interacting planets can have high-eccentricity or- 
bits that take them into previously undisturbed regions of 
the disc, is clearly the most challenging situation to model 
reliably with simple prescriptions. As shown, for example, 
in Figure [3] these phases can be important in shaping the 
long-term evolution of scattered systems, though such in- 
dividual examples do not mean that the statistical results 
obtained from a damping prescription might not be "good 
enough". Here, we have presented results from a direct hy- 
drodynamic simulation of a moderately large ensemble of 
systems, which we have used to extract statistical results 
for the eccentricity distribution and for the abundance of 
resonant planetary configuration. We find some differences 
between our hydrodynamic results and previously published 
results derived using approximate methods, though we de- 
fer for future work the question of whether other analytic 
prescriptions for planet-disc interactions would yield better 
agreement. 

Our main result is that scattering of random three- 
planet systems, from non-resonant initial conditions, 
yields significantly lower eccentricities (as compared to n- 
body runs) when the effects of hydrodynamics are included. 
This difference does not, however, arise because of any ma- 
jor change in the dynamics of scattering events that either 
eject planets, or result in planetary mergers. In these — 



dominant — channels, scattering in the hydrodynamic and 
n-body simulations leads to statistically very similar eccen- 
tricity distributions, differing (if at all) only in the fraction 
of highly eccentric planets. Rather, the lower eccentricities 
found in the hydrodynamic case occur because a modest 
number of scattering events evolve into three planet sys- 
tems that appear to be stable over an observationally in- 
teresting time scale (at least a Gyr). These systems, which 
tend to have lower than average eccentricities, are simply not 
present in runs without gas, whose end state after scattering 
is uniformly either one or two planet systems. 

Despite the lower eccentricities, we interpret our results 
as supporting the model of planetary scattering as the dom- 
inant physical mechanism giving rise to observed exoplanet 
eccentricity. Although gas does impact the outcome of scat- 
tering, it does not lead to large-scale circularization of plane- 
tary orbits. The differences between the hydrodynamic and 
gas- free simulations, in fact, are comparable to the varia- 
tions seen in scattering simulations that start with different 
planetary mass functions (Raymond et al. 2010). It would 



almost certainly be possible to find some set of plausible ini- 
tial conditions that would evolve hydrodynamically, under 
scattering, to match the observed distribution of massive 
planet eccentricity, given the freedom to fine tune the initial 
planetary spacing or disc dispersal time. Alternatively, and 
perhaps more naturally, the fraction of stable triple planet 
systems that are responsible for the lower eccentricity could 
likely be suppressed by starting with > 3. Given the 
wealth of such possibilities, our results do not constitute 
an argument against the physical importance of scattering, 
though they do imply that the relatively precise agreement 
between n-body calculations and observations is likely for- 
tuitous. 

The frequency of resonant systems of multiple planets 
is potentially a much more sensitive probe than the eccen- 
tricity distribution of the physics predating scattering. Un- 
der purely gravitational evolution, the occupation fraction 
of resonant systems mirrors the available phase space, and 
MMRs are uncommon (a few percent, 'Ray mond et al.|2008 ). 
In a laminar, viscous disc, conversely, almost arbitrarily high 
resonant fractions can be achieved ( Matsumura et al.|2010 ). 
Our results fall somewhere in between these extremes. The 
frequency of 2-planet resonances that we find is comparable 
to that found in pure n-body simulations. This is consistent 
with the conclusion, noted above, that the hydrodynamics 
and dynamics of scattering are very similar for those chan- 
nels that result in planetary ejection or merger. However, 
we also observe a high frequency of resonant chains of 3- 
planet systems. Out of 100 runs, 12 evolved into libration of 
the Laplace angle Lp l ■ This resonance trapping stabilized the 
systems against instability, pushing the number of unscat- 
tered systems in our runs to nearly 30%. Similar effects were 



seen in the n-body simulations of Raymond et al. (2010), 
where the dissipative role that is here played by gas was 
instead furnished by planet esimals. Observationally, some 
resonant systems are seen, but it is reasonably clear that 
their number is lower than theoretical expectations based 
on an extended period of evolution in a viscous disc. One 
possibility is that survival in resonance is frustrated at early 
times, when the disc is massive and turbulent fluctuations 



lead to stochastic torques on planets (Adams et al. 2008), 



such that resonance capture only becomes possible in a lim- 
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ited window just prior to disc dispersal. The quantitative 
details of such a scenario are, however, poorly known, and it 
is likely that our limited understanding of the mechanisms 
that can disrupt resonances is the dominant uncertainty in 
current models of early planetary system evolution. 

Finally we have explored the potential of future imag- 
ing to directly observe the disc structures formed during the 
epoch of planetary scattering, which include eccentric gaps 
and massive planets directly accreting from the disc gas. The 
predicted structures are distinctive, but the limited synchro- 
nization between the onset of scattering and the dispersal 
of the disc means that they should only rarely be present 
in a sample of transition discs. We estimate an abundance 
that might be of the order of a few percent. A blind search 
of transition discs is thus unlikely to uncover any systems 
with ongoing scattering, though in principle it may be possi- 
ble to exploit the large-scale asymmetries generated during 
scattering to identify possible candidates. 



ACKNOWLEDGMENTS 

Our thanks to the referee for a quick and constructive report. 
This work was performed in part using the Darwin Super- 
computer of the University of Cambridge High Performance 
Computing Service (http://www.hpc.cam.ac.uk/), provided 
by Dell Inc. using Strategic Research Infrastructure Fund- 
ing from the Higher Education Funding Council for Eng- 
land. PJA acknowledges support from NASA, under award 
NNX09AB90G and NNX11AE12G from the Origins of Solar 
Systems and Astrophysics Theory programs, and from the 
NSF under award AST-0807471. 



REFERENCES 

Adams F. C, Laughlin C, 2003, Icarus, 163, 290 
Adams F. C, Laughlin C, Bloch A. M., 2008, ApJ, 683, 
1117 

Aksnes K., 1988, in A. D. Roy ed.. Long-term Dynami- 
cal Behaviour of Natural and Artificial N-body Systems 
General formulas for three-body resonances., pp 125-139 

Beauge C, Giuppone C. A., Ferraz-Mello S., Michtchenko 
T. A., 2008, MNRAS, 385, 2151 

Butler R. P., Wright J. T., Marcy G. W., Fischer D. A., 
Vogt S. S., Tinney C. G., Jones H. R. A., Carter B. D., 
Johnson J. A., McCarthy C, Penny A. J., 2006, ApJ, 646, 
505 

Chambers J. E., 1999, MNRAS, 304, 793 
Chambers J. E., Wetherih G. W., Boss A. P., 1996, Icarus, 
119, 261 

Chatterjee S., Ford E. B., Matsumura S., Rasio F. A., 2008, 
ApJ, 686, 580 

D'Angelo G., Lubow S. H., Bate M. R., 2006, ApJ, 652, 
1698 

Fabrycky D. C, Murray-Clay R. A., 2010, ApJ, 710, 1408 

Fasano G., Franceschini A., 1987, MNRAS, 225, 155 

Ford E. B., Rasio F. A., 2008, ApJ, 686, 621 

Gladman B., 1993, Icarus, 106, 247 

Juric M., Tremaine S., 2008, ApJ, 686, 603 

Kley W., 1999, MNRAS, 303, 696 

Kley W., Peitz J., Bryden G., 2004, A&A, 414, 735 



Lee A. T., Thommes E. W., Rasio F. A., 2009, ApJ, 691, 
1684 

Lee M. H., Peale S. J., 2002, ApJ, 567, 596 
Lin D. N. C, Ida S., 1997, ApJ, 477, 781 
Lubow S. H., Ogilvie G. I., 2001, ApJ, 560, 997 
Lubow S. H., Pringle J. E., 2010, MNRAS, 402, L6 
Marchal C, Bozis G., 1982, Celestial Mechanics, 26, 311 
Marzari F., Baruteau C, Schoh H., 2010, A&A, 514, L4 
Marzari F., Nelson A. F., 2009, ApJ, 705, 1575 
Marzari F., Weidenschilling S. J., 2002, Icarus, 156, 570 
Masset F., 2000, A&AS, 141, 165 

Matsumura S., Thommes E. W., Chatterjee S., Rasio F. A., 

2010, ApJ, 714, 194 
Moeckel N., Raymond S. N., Armitage P. J., 2008, ApJ, 

688, 1361 

Moorhead A. V., Adams F. C, 2005, Icarus, 178, 517 
Murray CD., Dermott S. F., 1999, Solar system dynamics 
Owen J. E., Ercolano B., Clarke C. J., 2011, MNRAS, 412, 
13 

Owen J. E., Ercolano B., Clarke C. J., Alexander R. D., 
2010, MNRAS, 401, 1415 

Prince P. J., Dormand J. R., 1981, Journal of Computa- 
tional and Applied Mathematics, 7, 67 

Rasio F. A., Ford E. B., 1996, Science, 274, 954 

Raymond S. N., Armitage P. J., Gorelick N., 2010, ApJ, 
711, 772 

Raymond S. N., Barnes R., Armitage P. J., Gorelick N., 

2008, ApJ, 687, L107 
Rivera E. J., Laughhn G., Butler R. P., Vogt S. S., Haghigh- 

ipour N., Meschiari S., 2010, ApJ, 719, 890 
Snellgrove M. D., Papaloizou J. C. B., Nelson R. P., 2001, 

A&A, 374, 1092 
Terquem C, Papaloizou J. C. B., 2000, A&A, 360, 1031 
Thommes E. W., Lissauer J. J., 2003, ApJ, 597, 566 
Weidenschilling S. J., Marzari F., 1996, Nature, 384, 619 
Wolk S. J., Walter F. M., 1996, AJ, 111, 2066 



APPENDIX A: THE EFFECT OF 
ZERO-INCLINATION INITIAL CONDITIONS 

The FARGO and mercury scattering calculations reported 
in this paper employ exactly planar planetary initial condi- 
tions. This is done for reasons of computational economy: 
3-dimensional hydrodynamic simulations with the same in- 
plane resolution as the 2D runs would require around two 
orders of magnitude greater resources. For many purposes, 
moreover, assuming planarity is not a bad approximation. 
Prior n-body scattering calculations, run using the same 
mass distribution as ours, show that the median inclination 
of surviving planets is only ^ 10° (Raymond et al. 2010). 



This means that the initial phase of scattering typically in- 
volves inclinations i ^ h/r, or smaller, and accordingly we 
expect that any hydrodynamic interaction with the disc will 
be dominated by planar effects (eccentricity damping, accre- 
tion) rather than fully 3D ones (such as warp excitation). 
As noted in Section |3] however, the assumption of exact 
coplanarity does introduce a technical subtlety, since the fi- 
nal /(e) that results from scattering in planar systems is 
not the same as that which results from otherwise identical 
initial conditions with a small (of the order of one degree) 
non-zero inclination. Because of this subtlety, we compare 
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Figure Al. The cumulative eccentricity distributions, following 
scattering, obtained from purely n-body simulations of 3-planet 
systems initialized with different inclinations. The highest mean 
eccentricities result from initial conditions in which all three plan- 
ets had inclinations chosen randomly and uniformly in the interval 
i = [0, 1°] (blue curve). Lower eccentricities result from choosing 
i = [0, 0.1°] (green), or i = [0, 0.01°] (red). The black curve shows 
results for a strictly planar geometry. In all cases, only surviving 
planets from those systems that went unstable during the inte- 
gration are plotted. 



The results for the eccentricity distribution are shown 
in Figure |A1| We find that the ensembles initialized with 
i = [0, 0.1°] and i = [0, 0.01°] are similar to, but statistically 
distinct from, the zero inclination limit (K-S p-values 0.049 
and 0.023), whereas inclinations of the order of 1 degree are 
sufficient to result in markedly larger final eccentricities, dis- 
tinct from all other distributions. To an order of magnitude, 
the condition for scattering behavior that is distinct from 
the 3D limit in these three planet systems is thus that the 
initial inclinations should satisfy i < 0.1°. 

Could real planetary systems be flat enough that 
their scattering dynamics are effectively two-dimensional? 
It seems unlikely. Although a gas disc damps the inclination 
of a single planet on a short time scale, it will damp the or- 
bit until it is coincident with the local angular momentum 
vector of the gas, which will not be exactly constant with 
radius. A warp whose magnitude (over several AU) exceeds 
the threshold Ai = 0.1° seems almost inevitable (Terquem 



k Papaloizou 2000 Lubow k Pringle 2010). Even if the 



gas disc were perfectly fiat, any temporary resonances be- 
tween the planets would likely excite their mutual inclina- 
tion ( Thommes Lissauer|[2003 ) to the level where subse- 
quent scattering would occur in the three-dimensional limit. 
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our hydrodynamic simulations not to the same n-body runs 
that match the observed eccentricity distribution, but rather 
to their zero inclination analogues. 

We emphasize that the point of this paper is not to ar- 
gue that the result of the planar calculation is physically 
realistic. Rather, from the fact that in 2D the effect of hy- 
drodynamics on the final eccentricity distribution is small, 
we infer that in 3D a hydrodynamic calculation would yield 
results similar to the observed distribution. That said, a 
physical effect is in principle possible if the gas disc were 



able to damp mutual inclinations ( [Lubow Og ilvie 2001 



Marzari & Nelson 2009) to the point where the final sys- 



tem was effectively planar. To quantify this possibility, we 
have run a set of n-body mercury experiments to identify 
how small the initial inclination must be for the system to 
scatter as if it were exactly planar. The setup is identical 
to that used earlier. We start with three planets on circu- 
lar orbits, separated by K = 4.0 mutual Hill radii (Equa- 
tion[T]), with the inner planet at 3 AU. The mass spectrum is 
f{m) oc m~^'^ between 0.3 Mj and 5 Mj, and the planetary 
radius is 1.3 Rj. We run four ensembles of 500 simulations, 
with the planetary inclinations chosen randomly and uni- 
formly in the range i = [0,1°], i = [0,0.1°], i = [0,0.01°] 
and i = 0. The systems are integrated for 15 Myr, after 
which we plot the eccentricity distribution of the surviving 
planets in the unstable systems (at this epoch, this is about 
80% of the total). As before, we discard double systems that 
are not Hill stable, which removes only ~ 10 runs from each 
set. 
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